//============================================================================
// Name : monte.cpp
// Author : Nazmus Saquib
// Version :
// Copyright :
// Description : BCA code in C++
//============================================================================
#include “Globals.h”
#include “scoef.h”
#include “rstop.h”
#include “muscle_bca.h”
#include “vectorGeometry.h”
const int numXBin = 100;
const int numIons = 1;
int xBin[numXBin+1] = {0};
int selectedXBin = 0;
double rho[4]={0.0};
int n[]={0,0,0,0};
double e0kev, m1, cw, ed, latticeConst;
int z1, hn, iy, nowout;
double dx[3] = {0.0};
int zt[3][7] = {0};
double mt[3][7] = {0.};
double t[3][7] = {0.};
int Layer; //use this for number of layers rather than the L in the program;
//remember this cannot be set to 0. Other loops/arrays may start from
//0, but keep this >= 1;
//the arrays..
double my[3][7]= {0},
ai[3][7]= {0},
fi[3][7]= {0},
ec[3][7]= {0},
io[3][7]= {0},
k[3]= {0},
kl[3][7]= {0};
double vf[3][7]= {0},
mu[3]= {0},
ioniz[3]= {0},
h[3]= {0},
xx[3]= {0},
m2[3]= {0.0},
z2[3]= {0},
c[3]= {0},
epsbk[3]= {0},
arho[3]= {0};
double a[3]= {0},
f[3]= {0},
lm[3]= {0},
pmax[3]= {0},
fd[3]= {0},
kd[3]= {0},
sbk[3]= {0},
lf[3]= {0},
yy[8]= {0};
double se[3][1000]= {0.0},
seo[1000]= {0.0},
epsdg[3]= {0.0};
double ls = 0., lo = 0., maximum = 0.;
double xsum = 0, x2sum = 0, x3sum = 0, x4sum = 0, plsum = 0, pl2sum = 0;
double avex = 0, vari = 0, sigma = 0, v = 0, v2 = 0, Gamma = 0, beta = 0, y = 0, avepl
= 0, sigpl = 0, avecol = 0;
int i = 0, j = 0;
int ib = 0, it = 0;
double eb = 0.0, et = 0.0;
int icsum = 0;
double y2sum = 0, xy2sum = 0, x2y2su = 0, y4sum = 0;
double tau = 0;
int iii = -1;
//my dummy vars
int w; double q;
//my customized data structures or vars
scoefData scoefz1;
rstopData rstp;
//these vars (continuing from this line) are added later as needed, vars that were not
declared in the initialization of trim85 but showed up in the middle.
double e0, ef;
int iz;
double alfa, alpha;
double tmin, da;
double L0;
int iz1;
double ee;
int izt;
double nh;
double epso;
int ih;
double e;
//boolean for keeping track of a particle being transmitted or backscattered
int transmitted = 0;
int backscattered = 0;
//boolean for keeping track of channeling
int insideChannel = 1;
//boolean for determining whether to scatter from neighbor atoms
int neighborFlag = 0;
//main function
int main()
{
cout << “Initializing” << endl;
Initialize(“IronInput.dat”);
calculateAvgMassOfLayer();
getStoppingForTarget();
setInitialConditions();
MonteCarlo();
/* *********************** Testing Ground for arrays and variables
******************** */
for (w = 1; w <= 8; w++) cout << “\t “ << yy[w];
for (w = 1; w <= 3; w++) cout << “\n\t “ << xx[w];
cout << “CW or L0: “ << L0;
//(trouble! The values of n[] elements are changing like crazy)
cout << endl << n[0] << “\t” << n[1] << “\t” << n[2] << “\t” << n[3] << endl;
cout << endl << rho[0] << “\t” << rho[1] << “\t” << rho[2] << “\t” << rho[3] <<
“\t” << rho[4] << endl;
//test z2[], m2[] arrays
for (w = 1; w <= Layer; w++) cout << “\t “ << z2[w];
for (w = 1; w <= Layer; w++) cout << “\n\t “ << m2[w];
//test the se values
cout << endl << se[1][1000] << “ “ << se[2][1000] << “ “ << se[3][1000] << endl;
//test the se[] values by printing 10th element from the array (to be implemented)
//for (w = 1; w <= 100; w++)
//{
//for (int sss = 1; sss <= 50; sss++)
//{
//cout << “\t” << mpart[w][sss];
//}
//}
/* *********************** End of Testing Ground ************************** */
cout << endl;
//Now print values of each bin in xBin
//for (w = 1; w <= numXBin; w++) cout << xBin[w] << “,”;
//print final x values
cout << endl;
cout << “Number of Backscattered Ions: “ << ib << endl;
cout << “Number of Transmitted Ions: “ << it << endl;
return 0;
}
//end of main
//Helper functions follow from here
//initialize variables and setup each layer over here
void Initialize(char* filename)
{
//Read command file
ifstream instream;
instream.open(filename,ios::in);
cout << filename;
if (!instream)
{
exitOnError(filename);
}
instream >> e0kev >> z1 >> m1 >> latticeConst >> hn >> cw >> ed >> iy >> nowout;
instream >> dx[1] >> rho[1];
instream >> zt[1][1] >> mt[1][1] >> t[1][1];
instream >> n[1];
instream >> dx[2] >> rho[2];
instream >> zt[2][1] >> mt[2][1] >> t[2][1];
instream >> n[2];
instream >> dx[3] >> rho[3];
instream >> zt[3][1] >> mt[3][1] >> t[3][1];
instream >> n[3];
for (w = 0; w <= 3; w++)
{ if (n[w] != 1) n[w] = 1; }
instream.close();
Layer = 3;
//done with primary (crude) setup, off to reading scoef.dat file
iz = z1;
getstructScoef(iz, scoefz1, “scoef1.dat”); //used to get pcoef data yy in trim85
//so we put the values in yy here immediately
for (w = 1; w <= 8; w++)
{ yy[w] = scoefz1.pcoef[w]; }
//note: although yy turns out to be just a dummy array
e0 = e0kev*1000; //convert to ev.
if (ed == 0.) ed = 25.0;
ef = Max(5.0, e0kev*0.1);
//alfa = angle of incidence, alpha = radian of alfa
alfa = 0.;
alpha = alfa*Pi/180;
tmin = 5.0;
tau = 0.0;
da = 3.0;
if (iy == 0) iy = 16381;
//now calculate the total depth of each layer = xx(l), and grid spacing cw
xx[1] = dx[1];
for (w = 2; w <= 3; w++) { xx[w] = dx[w] + xx[w-1]; }
if (cw == 0) cw = 0.01*xx[3];
L0 = cw;
//take care of any custom variable or struct I made
rstp.vfermi = 0.0; //set this to 0. for the time being (see log for explanation)
//now, off to avg mass of layer in the next void function
}
//avg mass and atomic number of each layer
void calculateAvgMassOfLayer()
{
for (int LL = 1; LL <= Layer; LL++)
{
int ii = n[LL];
for (w = 1; w <= ii; w++)
{
h[LL] = h[LL] + t[LL][w];
//cout << “testing here..” << rstp.vfermi << endl;
}
}
for (int LL = 1; LL <= Layer; LL++)
{
int ii = n[LL];
for (w = 1; w <= ii; w++)
{
t[LL][w] = t[LL][w]/h[LL];
m2[LL] = m2[LL] + t[LL][w] * mt[LL][w];
z2[LL] = z2[LL] + t[LL][w] * zt[LL][w];
//note: z2 wont be an integer once t has a value other than 1.0??
//so are we calculating average atomic number here? why?
}
}
//done with this, off to finding electronic stopping powers in the next function
}
void getStoppingForTarget()
{
iz1 = z1; ee = 0;
for (int LL = 1; LL <= Layer; LL++)
{
arho[LL] = rho[LL] * 0.6022/(m2[LL]);
mu[LL] = m1/(m2[LL]);
int ii = n[LL];
for (int nn = 1; nn <= ii; nn++)
{
//set rstp.se[1..1000] = 0. Clear it for the next loop
for (w = 1; w <= 1000; w++)
{ rstp.se[w] = 0.; }
izt = zt[LL][nn];
//calling getrstop now. units, lfctr, vfermi = 1 (doesnt matter)
getrStop(iz1, izt, e0kev, 1, 1., 1., rstp); //rstp defined in header
//set this anyway, though I dont calculate this in RSTOP
vf[LL][nn] = rstp.vfermi; //which is just set 0 in the struct def.
for (w = 1; w <= 1000; w++)
{
se[LL][w] = se[LL][w] + rstp.se[w] * t[LL][nn] * arho[LL];
}
}
}
//now, off to setting up initial conditions next function
}
void setInitialConditions()
{
nh = hn; //number of histories
for (int LL = 1; LL <= Layer; LL++)
{
a[LL] = 0.5292 * 0.8853 / ( pow(z1, 0.23) + pow(z2[LL], 0.23) );
//now calculate the mean flight path with the conditions given in trim85
f[LL] = a[LL] * m2[LL] / ( z1 * z2[LL] * 14.4 * (m1 + m2[LL] ) );
epso = e0 * f[LL];
epsdg[LL] = tmin * f[LL] * pow( (1.0 + mu[LL]) , 2) / (4.0 * mu[LL]);
fd[LL] = 0.01 * pow( z2[LL], (-7.0/3.0) );
kd[LL] = 0.1334 * pow ( z2[LL], (2.0/3.0) ) / sqrt( m2[LL] );
}
for (int LL = 1; LL <= Layer; LL++)
{
int ii = n[LL];
for (w = 1; w <= ii; w++)
{
my[LL][w] = m1/mt[LL][w];
ec[LL][w] = 4.0 * my[LL][w] / pow( ( 1.0 + my[LL][w] ), 2);
ai[LL][w] = 0.5292 * 0.8853 / ( pow(z1,0.23) + pow(zt[LL][w],0.23) );
fi[LL][w] = ai[LL][w] * mt[LL][w] / ( z1 * zt[LL][w] * 14.4 * ( m1 + mt[LL][w] ) );
}
}
cout << “Setup finished. Starting Monte Carlo Loops..”;
//off to monte carlo loop in next routine
}
void MonteCarlo()
{
//custom variables and arrays for this section
double e;
double cosin = 0.0, siny = 0.0, sine = 0.0, cosy = 0.0;
double pl = 0.0;
int ic;
int LL = 1;
double eps;
double eeg;
double p;
double b;
int ie, ia; //not sure if ie or ia should be double. They are used to access elements
of the m[][] array at some point
double see;
double dee;
double s2, c2, ct, st;
double r, rr;
double ex1, ex2, ex3, ex4;
double v, v1;
double fr, fr1;
double q;
double roc, sqe;
double cc, aa, ff;
double delta, co;
double den;
double phi, psi;
double x1;
int ip;
//variables for crystal calculations
double sep = 0.0;
int ionCounter = 0;
double p1, p2;
double Theta = 0.0;
double rTheta = 0.0, rPhi = 0.0;
double thetaThreshold = 0.5;
double crystalMuonY = 0.0, crystalMuonZ = 0.0;
//amount of translations in y and z axes
double translationY = 0.0, translationZ = 0.0;
//Scatter Plot variables
const int numScatterPlotBins = 5;
double scatterPlot[numScatterPlotBins + 1] = {0.0};
//Vector declarations
//lattice constant of Target (input from file)
double latticeConstant = latticeConst;
Vector3 ions[4];
Vector3 ionsOrdered[4];
Vector3 neighborIonsBCC[14];
Vector3 neighborIonsFCC[14];
Vector3 ionTranslationY(0,translationY,0);
Vector3 ionTranslationZ(0,0,translationZ);
Vector3 unitzplus, unitzminus, unityplus, unityminus, unitxplus, unitxminus;
unitzplus.z = 1;
unitzminus.z = -1;
unityplus.y = 1;
unityminus.y = -1;
unitxplus.x = 1;
unitxminus.x = -1;
Vector3 initialDirection(1,0,0);
Vector3 d;
d.x = 1;
d *= latticeConstant/2;
Vector3 lx(latticeConstant/2, 0, 0);
Vector3 ly(0, latticeConstant/2, 0);
Vector3 lz(0, 0, latticeConstant/2);
Vector3 lambda;
Vector3 lambdaPrime;
Vector3 pVector;
Vector3 pUnitVector;
Vector3 sepVector;
Vector3 Di;
Vector3 DiPrev;
//Vector3 DiPrevToDi;
Vector3 delX;
Vector3 delX1, delX2;
Vector3 dummy1, dummy2, temp;
Vector3 scatterIonPos;
//initialize the random number generator
srand(time(NULL));
//open file to write output
ofstream outStream;
outStream.open(“coords0.txt”);
if(outStream.fail())
{
exitOnError(“Could not open Output file”);
}
//write basic information in the output file
//number of ions, ion energy, total depth, depth of each layer
outStream << nh << “\t” << e0kev << “\t” << xx[Layer] << “\t” << dx[1] << “\t”
<< dx[2] << “\t” << dx[3] << endl;
//open scatter plot file to write current Y and Z coordinates of muons at designated
intervals
ofstream scatterStream;
scatterStream.open(“scatterOut1.txt”);
if(scatterStream.fail())
{
exitOnError(“Could not open Scatter Plot Output file”);
}
//open range distribution file to write final X coordinates of muons
ofstream rangeStream;
rangeStream.open(“rangeOut1.txt”);
if(rangeStream.fail())
{
exitOnError(“Could not open Range Distribution Output file”);
}
//Open general information dump file
ofstream infoStream;
infoStream.open(“info1.txt”);
if(infoStream.fail())
{
exitOnError(“Could not open general information Output file”);
}
//Entering the target
//First set up for the top layer
for (ih = 1; ih <= nh; ih++)
{
avex = xsum / Max(1.0, (float)(ih - ib - it - 1));
if (talk == 2) cout << “Average ion range so far: “ << avex << “ angstroms.”
<<endl;
if (talk > 2) cout << “Now starting ion number “ << ih << endl;
e = e0;
//set scatterPlot array (the intervals)
scatterPlot[0] = 10;
scatterPlot[1] = 30;
scatterPlot[2] = 50;
scatterPlot[3] = 100;
scatterPlot[4] = 150;
scatterPlot[5] = 200;
/*
for(int ccc = 0; ccc <= numScatterPlotBins; ccc++)
{
cout << “scatter plot “ << ccc << “: “ << scatterPlot[ccc] << endl;
}
*/
//Initial Ion Positions
ions[0] = d + ionTranslationZ + unitzminus * (latticeConstant/2);
ions[1] = d + ionTranslationZ + unitzplus * (latticeConstant/2);
ions[2] = d * 2 + ionTranslationY + unityminus * (latticeConstant/2);
ions[3] = d * 2 + ionTranslationY + unityplus * (latticeConstant/2);
//Create a polygon that resides on the lateral axes.
//The points are put on anticlockwise order, which is important for
//testing whether the test point lies on this polygon
ionsOrdered[0] = ions[0];
ionsOrdered[1] = ions[2];
ionsOrdered[2] = ions[1];
ionsOrdered[3] = ions[3];
//generate random theta and phi angles.
rTheta = ( (double)rand()/((double)(RAND_MAX)+(double)(1)) ) * thetaThreshold;
rPhi = ( (double)rand()/((double)(RAND_MAX)+(double)(1)) ) * 2 * Pi;
//find corresponding x, y and z components of the direction vector.
//radius of the direction vector is 1
initialDirection.x = cos(rTheta);
initialDirection.z = sin(rTheta) * cos(rPhi);
initialDirection.y = sin(rTheta) * sin(rPhi);
//set the counter to 0 for a new ion
ionCounter = 0;
pl = 0.0;
ic = 0;
//set initial DiPrev - the origin
DiPrev.x = 0.0; DiPrev.y = 0.0; DiPrev.z = 0.0;
//set initial lambda, the direction of motion. Normalize it.
lambda.clear();
lambda = initialDirection;
lambda.normalize();
//set initial delX to origin
delX.clear();
//clear the dummy delX vectors, set them to origin
delX1.clear();
delX2.clear();
dummy1.clear();
dummy2.clear();
temp.clear();
LL = 1;
//set transmitted and backscattered to false
transmitted = 0;
backscattered = 0;
//set channeling to true
insideChannel = 1;
neighborFlag = 0;
//write the initial coordinates to the output file
outStream << DiPrev.x << “\t” << DiPrev.y << “\t” << DiPrev.z << endl;
//cout << endl << “Initial: “;
//cout << DiPrev.x << “\t” << DiPrev.y << “\t” << DiPrev.z << endl;
//cout << “r1 = “ << r1 << endl;
//cycle for each collision until the energy of the particle becomes too low, or
the particle backscatters, or it goes out of the last layer (transmission)
//needs a do while loop here,
//which I will mention as the ‘mother loop’ from now.
do
{
ic = ic + 1;
eps = e * f[LL];
eeg = sqrt(eps*epsdg[LL]);
//pmax[LL] = a[LL] / (eeg + sqrt(eeg) + 0.125 * pow( eeg, 0.1) );
pmax[LL] = sqrt(3) * (latticeConstant / 2) * 0.7;
//Calculate impact parameter and choose the atom to scatter from.
//Do this for ion pairs 0,1 and 2,3.
if (ionCounter == 0)
{
delX1 = ions[0] - DiPrev;
delX2 = ions[1] - DiPrev;
dummy1 = delX1 % lambda;
dummy2 = delX2 % lambda;
p1 = sqrt( dummy1.scalarProduct( delX1 % lambda ) );
p2 = sqrt( dummy2.scalarProduct( delX2 % lambda ) );
if(p2 > p1)
{ //swap ion ordering
temp = ions[0];
ions[0] = ions[1];
ions[1] = temp;
//cout << “Vertical Ions swapped” << endl;
}
}
if (ionCounter == 2)
{
delX1 = ions[2] - DiPrev;
delX2 = ions[3] - DiPrev;
dummy1 = delX1 % lambda;
dummy2 = delX2 % lambda;
p1 = sqrt( dummy1.scalarProduct( delX1 % lambda ) );
p2 = sqrt( dummy2.scalarProduct( delX2 % lambda ) );
if(p2 > p1)
{ //swap ion ordering
temp = ions[2];
ions[2] = ions[3];
ions[3] = temp;
//cout << “Horizontal Ions swapped” << endl;
}
}
//now calculate impact parameter
if(neighborFlag == 0)
{
delX = ions[ionCounter] - DiPrev;
dummy1 = delX % lambda;
p = sqrt( dummy1.scalarProduct( delX % lambda ) );
//find impact parameter vector and it’s unit vector
pVector = dummy1 % lambda;
pUnitVector = pVector.unit();
}
else if(neighborFlag == 1)
{
double impact[13] = {0.0};
double radial[13] = {0.0};
double S[13] = {0.0};
double sumS = 0.0;
double Probability[13] = {0.0};
double rnd_candidate = 0.0;
int selected_candidate = -1;
for(int ncount = 0; ncount <= 13; ncount++)
{
delX = neighborIonsBCC[ncount] - DiPrev;
dummy1 = delX % lambda;
p = sqrt( dummy1.scalarProduct( delX % lambda ) );
//find impact parameter vector and it’s unit vector
pVector = dummy1 % lambda;
pUnitVector = pVector.unit();
//scatterIonPos = neighborIonsBCC[ncount];
impact[ncount] = p;
radial[ncount] = delX.magnitude();
S[ncount] = 1 / ( pow(impact[ncount],2) * radial[ncount] );
sumS += S[ncount];
//general scheme of selecting the neighbor ion
// if(p < pmax[LL])
// {
// scatterIonPos = neighborIonsBCC[ncount];
// //cout << “Neighbor “ << ncount << “ is selected” << endl;
// break;
// }
}
for(int ncount = 0; ncount <= 13; ncount++)
{
Probability[ncount] = S[ncount] / sumS;
}
//print out the probability array
cout << endl;
for(int ncount = 0; ncount <= 13; ncount++)
{
//cout << Probability[ncount] << “ “;
infoStream << Probability[ncount] << “ “;
}
infoStream << endl;
//cout << endl;
//random number between 0 and 1
rnd_candidate = ( (double)rand()/((double)(RAND_MAX)+(double)(1)) );
//cout << “rand_candidate: “ << rnd_candidate << endl;
//choose the candidate for scattering
selected_candidate = whichNonUniformBin(rnd_candidate, Probability, 13);
//cout << “Ion: “ << ih << “\tSelected Candidate: “ << selected_candidate <<
endl;
infoStream << “Ion: “ << ih << “\tSelected Candidate: “ << selected_candidate
<< endl;
//assign the scatterIonPos variable to the selected neighbor
scatterIonPos = neighborIonsBCC[selected_candidate];
//cout << “Scattering Ion Position: “; scatterIonPos.printVector();
//find the essential quantities for the selected neighbor
delX = neighborIonsBCC[selected_candidate] - DiPrev;
dummy1 = delX % lambda;
p = sqrt( dummy1.scalarProduct( delX % lambda ) );
//find impact parameter vector and it’s unit vector
pVector = dummy1 % lambda;
pUnitVector = pVector.unit();
}
//find eps and b using fi[LL][nn], using nn that I was supposed to find from above
//here im deliberately using nn = 1
eps = fi[LL][1] * e;
b = p / ai[LL][1];
if (eps > 10) //rutherford scattering
{
s2 = 1.0 / (1.0 + (1.0 + b * (1.0 + b)) * pow((2.0 * eps * b), 2) );
c2 = 1.0 - s2;
ct = 2.0 * c2 - 1.0;
st = sqrt(1.0 - ct * ct);
}
else //magic formula
{
r = b;
rr = -2.7 * log(eps * b);
if (rr >= b)//note >= sign instead < in trim85
{
rr = -2.7 * log(eps * rr);
if (rr >= b)//note >= sign instead < in trim85
{
r = rr;
}
}
//do while loop that replaces line 330 loop
do
{
ex1 = 0.18175 * exp(-3.1998 * r);
ex2 = 0.50986 * exp(-0.94229 * r);
ex3 = 0.28022 * exp(-0.4029 * r);
ex4 = 0.028171 * exp(-0.20162 * r);
v = (ex1 + ex2 + ex3 + ex4) / r;
v1 = -(v + 3.1998 * ex1 + 0.94229 * ex2 + 0.4029 * ex3 + 0.20162 * ex4) / r;
fr = b * b / r + v * r / eps - r;
fr1 = -b * b / (r * r) + (v + v1 * r) / eps - 1.0;
q = fr / fr1;
r = r - q;
}
while( (Abs(q / r)) > 0.001 );
roc = -2.0 * (eps - v) / v1;
sqe = sqrt(eps);
//5 parameter magic scattering calculation
//below is for universal potential
cc = (0.011615 + sqe) / (0.0071222 + sqe);
aa = 2.0 * eps * (1.0 + (0.99229 / sqe) ) * ( pow(b, cc) );
ff = ( sqrt(aa * aa + 1.0) - aa) * ( (9.3066 + eps) / (14.813 + eps) );
delta = (r - b) * aa * ff / (ff + 1.0);
co = (b + delta + roc) / (r + roc);
c2 = co * co;
s2 = 1.0 - c2;
ct = 2.0 * c2 - 1.0;
st = sqrt(1.0 - ct * ct);
}
Theta = acos(ct);
//we are done finding theta (in CM system). So calculate all other quantities.
//find separation and the separation vector.
phi = (Pi - Theta) / 2;
sep = p / tan(phi);
sepVector = lambda * sep;
//find theta in laboratory frame - psi
psi = atan(st / (ct + my[LL][1] ) );
//note: change my[LL][1] to my[LL][nn] when the above section is fixed.
if (psi < 0 ) psi = psi + Pi; //should I do this for crystals?
//find Di, the scattering point vector
Di = DiPrev + delX + pVector - sepVector;
//find new direction of motion
lambdaPrime = lambda * cos(psi) + pUnitVector * sin(psi);
lambdaPrime.normalize();
//find length of step, ls = distance of Di from DiPrev
ls = Di.getDistance(DiPrev);
//find energy lost due to electronic stopping, dee
ie = (int)(e/e0kev+0.5); //should it be 0.5? or less so that ie <=1000?
see = se[LL][ie];
if (e < e0kev) see = se[LL][1] * sqrt(e/e0kev);
dee = ls * see;
// den = energy transferred to recoil
den = ec[LL][1] * s2 * e; //note: I am using ec[LL][1] here instead of [LL][nn].
//cout << “den = “ << den << “, dee = “ << dee << endl;
infoStream << “den = “ << den << “, dee = “ << dee << endl;
e = e - den - dee;
//cout << endl<< “current ion energy: “ << e << endl;
if (dee > maximum) maximum = dee;
pl = pl + ls - tau;
//write the ion position to output file
outStream << Di.x << “\t” << Di.y << “\t” << Di.z << endl;
if((ic%30)==0)
{
//cout << Di.x << “\t” << Di.y << “\t” << Di.z << endl;
}
//determine which scatter plot Di’s x value belongs to.
//output the y and z coordinates to scatter plot file accordingly.
for(int cc = 0; cc <= numScatterPlotBins; cc++)
{
if( Di.x >= (scatterPlot[cc]) )
{
crystalMuonY = Di.y - translationY;
crystalMuonZ = Di.z - translationZ;
//cout << “scatter plot “ << ct << “: “ << scatterPlot[ct] << endl;
scatterStream << scatterPlot[cc] << “\t” << crystalMuonY << “\t” <<
crystalMuonZ << endl;
//set scatterPlot[ss] to a big number
scatterPlot[cc] = 10000;
//break out of this for loop
break;
}
}
//determine if Di is in the channeling region.
//insideChannel = Di.isInsidePolygon(ionsOrdered, 4);
//break out of parent loop if not inside the channel
if(insideChannel == 0)
{
cout << “Ion “ << ih << “ is out of Channel” << endl;
cout << Di.x << “\t” << Di.y << “\t” << Di.z << endl;
break;
}
//determine which layer the next collision will be in
if (Di.x < 0.0) //particle is backscattered
{
backscattered = 1; //cout << “ion number “ << ih << “ backscattered.” << endl;
infoStream << “ion number “ << ih << “ backscattered.” << endl;
ib = ib + 1;
eb = eb + e;
break; //break out of ‘mother do loop’ and continue with the next session of
for loop.
}
//here we set the current or next layer the particle will be in
for (w = 1; w <= Layer; w++)
{
if ( (Di.x <= xx[w]) && (w == 1) )
{
LL = 1;
//cout << endl<<”ion is in layer “ << LL << endl;
break;
//break out of this For loop and go check if the particle is transmitted.
}
else if ( (Di.x <= xx[w]) && (Di.x > xx[w-1]) )
{
LL = w;
//cout << “ion is in layer “ << LL <<endl;
break; //break out of this For loop and go check if the particle is
transmitted.
}
}
//now, check for particle transmission, i.e. whether the particle went out of the
last layer.
if(Di.x >= xx[Layer])
{
//particle is transmitted, take care of appropriate variables and break
transmitted = 1;//cout << “ion number “ << ih << “ transmitted.” << endl;
it = it + 1;
et = et + e;
ia = 57.295779 * acos(cosin) / da + 1.0;
ie = 100 * e / e0 + 1.0;
//m[ie][ia] = m[ie][ia] + 1;//note: how is this possible? ie and ia should be
integers in order to access the elements of the array m[][]. But we calculate them as
doubles here!
break; //break out of the ‘mother’ do loop
}
//now take care of ionCounter and other variables for the next scattering
if(neighborFlag == 0)
{
if(ionCounter == 3)
{
//ionCounter = 0;
neighborFlag = 1;
scatterIonPos = ions[3];
//cout << endl << “got out of first two layers” << endl;
}
else
{
ionCounter++;
}
}
if(neighborFlag == 1)
{
//cout << “Updating Neighbor Ions” << endl;
neighborIonsBCC[0] = scatterIonPos + lx - ly + lz;
neighborIonsBCC[1] = scatterIonPos + lx + ly + lz;
neighborIonsBCC[2] = scatterIonPos + lx + ly - lz;
neighborIonsBCC[3] = scatterIonPos + lx - ly - lz;
neighborIonsBCC[4] = scatterIonPos + lx * 2;
neighborIonsBCC[5] = scatterIonPos - lx - ly - lz;
neighborIonsBCC[6] = scatterIonPos - lx + ly - lz;
neighborIonsBCC[7] = scatterIonPos - lx + ly + lz;
neighborIonsBCC[8] = scatterIonPos - lx - ly + lz;
neighborIonsBCC[9] = scatterIonPos + lz * 2;
neighborIonsBCC[10] = scatterIonPos - lz * 2;
neighborIonsBCC[11] = scatterIonPos + ly * 2;
neighborIonsBCC[12] = scatterIonPos - ly * 2;
neighborIonsBCC[13] = scatterIonPos - lx * 2;
}
//set DiPrev to Di
DiPrev = Di;
//update current lambda to lambdaPrime
lambda = lambdaPrime;
//now the while condition of the mother do loop checks if the particle has lesser
energy than our lowest energy limit, ef.
}
while(e > ef);
//since we are out of the mother do loop now, the particle must have come to a
stop. So, increase the final particle distributions if the particle has not been
transmitted or backscattered.
if( ((transmitted == 0)) && ((backscattered == 0)) && ((insideChannel == 1)) )
{
ip = (int)(pl/cw + 1.0);
if(ip > 100) ip = 100;
//ipl[ip] = ipl[ip] + 1;
xsum = xsum + Di.x;
//my own bin function
selectedXBin = whichBin(Di.x, xx[Layer], numXBin);
//write the selected bin to the output file
//outStream << selectedXBin << endl;
//cout << x << endl << xx[Layer] << endl << numXBin << endl << selectedXBin;
xBin[selectedXBin] = xBin[selectedXBin] + 1;
//print final x value for plotting histogram
//cout << “ion “ << ih << “ final x: “ << Di.x << endl;
infoStream << “ion “ << ih << “ final x: “ << Di.x << endl;
//cout << Di.x << “,”;
rangeStream << Di.x << endl;
plsum = plsum + pl;
icsum = icsum + ic;
//ipl is the ion path length - the total avg. distance the ion travels
regardless of direction before it comes to stop
}
//that brings us to the end of one ion’s journey, now go to next ion by going back
to the for loop’s beginning..
}
//and this ends the monte carlo loop function. Take care of necessary structures and
variables that need to be cleared/deleted
outStream.close();
scatterStream.close();
infoStream << “Number of Backscattered Ions: “ << ib << endl;
infoStream.close();
rangeStream.close();
}
int whichNonUniformBin(double e, double arr[], int numBins)
{
double low = 0, high = arr[0];
if((e>=low) && (e<high))
return 0;
else
{
for(int i = 0; i < numBins; i++)
{
low += arr[i];
high += arr[i+1];
//cout << “\tlow: “ << low << “, high:” << high << endl;
if((e>=low) && (e<high))
return i;
}
}
}